# 进行逻辑回归学习 2023年3月3日16:31:37 
#*************         NHANES         *************
#*************      Analyze Code      *************
#************* Logistic 回归模型讲解  *************
#*************                        *************
#### 0 准备好环境####
library(tidyverse)
library(gtsummary)
library(haven)
# drop_na 函数，用于快速去掉 NA
library(tidyr)
library(survey)
setwd('G:/BaiduNetdiskDownload/NHANES')

#### 1. 常规Logistics模型 ####
##### 1.1 GLM 函数 #####
# Logit模型结果(是/否)所包含的信息要比定量结果少得多。因此，拟合一个合理的 Logit模型需要的数据比一个相似大小的线性模型要多。
#一个 Logit模型通常至少需要100个观测数据才能符合截距，另外每加预测因子（核心变量+协变量）还需要8-15个观测数据
# Logit 模型中Y变量的“是”或者“否”的占比也会影响到所需要的数据量，
# 通常可以以这个方式进行估计：100+10 * (#X)/(较小占比的Y的结果）
# e.g. #X = 5；Y-1-10%； 100+ 10*5/0.1 = 100+500 数量

# 生成 weight 变量，100 个男性 & 100 个女性
n = 100
weight.male <- rnorm(n, 75, 10)
weight.female <- rnorm(n, 55, 5)
weight <- c(weight.male, weight.female)

# 生成性别变量，100 个男性 & 100 个女性
sex.male <- factor(sample(c('male'), n, TRUE))
sex.female <- factor(sample(c('female'), n, TRUE))
sex <- c(sex.male, sex.female)


# 如果要随机生成 100 个男性和女性
# sex <- factor(sample(c('male', 'female'), n, TRUE))

# 生成身高变量
height <- weight + rnorm(n * 2, 0, 5) # 身高变量为 weight + 一个随机的、有点大的扰动项

# 构建简单的Logistic回归模型
# 注意：我们使用性别变量作为 Y，看下 weight、height 能否预测性别
# 要先把 Y 变为 0、1 的值, 更直观，不变的话也不会报错，会自动识别
sex.0.1 <- ifelse(sex == 'male', 1, 0)
data <- data.frame(height = height, weight = weight, sex = sex, sex.0.1 = sex.0.1)


m1 <- glm(
  sex.0.1 ~ weight + height, #
  data = data,
  family = binomial(link = "logit")
  # 这里也可以直接写 family = "bionomial", 因为 logit 链接是二项式家族的默认方法
)

# 系数输出提供了估计系数及其标准误差，再加上 Wald Z 统计量，即估计系数除以其标准误差。将其与标准的 Normal 分布进行比较，以获得在 Pr (> | z |)列中总结的双尾 p 值。
summary(m1)


# prz 越低 说明越有联系
coefficients(m1) 
exp(coefficients(m1)) # weight 前的系数表明，odds 的笔直随着weight每增加1，升高1.42 ，身高越大，Y取1（男性）的概率越大
confint(m1) #与线性回归一样，我们可以得到截距系数和斜率系数的95% 置信区间(为了得到其他水平，改变限制中的水平参数)。
fitted(m1)
AIC(m1) # AIC 越低代表模型的预测效果越好
predict(m1)

##### 1.2 使用 tbl_regression 优化输出 #####

m1_tbl_1 <-tbl_regression(m1)
m1_tbl_1


##### 1.3 可视化结果  #####
# ggplot(data, aes(x = weight, y = sex.0.1)) +
#   geom_jitter(height = 0.05) +
#   geom_line(aes(x = weight, y = fitted(m1)), 
#             col = "blue") +
#   labs(title = "Logistic Regression from Model")


